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Abstract 

We  explore  experimental  procedure?  for  comparing  the  capal, !li ties  of  complex  discrete  <  vent 
service  systems.  Instead  of  measuring  system  capability  by  analyzing  or  simulating  the  system 
with  a  constant  rate  of  arriving  work,  system  capability  is  measured  as  the  maximum  rate  of 
work  arrival  for  which  the  system  has  a  steady  state.  Hence,  we  seek  the  arrival  rate  which 
causes  the  system  to  be  at  full  capacity.  This  rate  is  arguably  the  best  indication  of  the  service 
system's  capability.    We  treat  both  work-conserving  and  non-work-conserving  service  systems, 

U-MUg,   iiauiuuiidi    ctnu   spccidaiicu   mcd-MUcs  ui  s^mciii    jit  1  lul  liidlit-c. 

Introduction 

As  industrial  engineers,  applied  probabilists,  simulationists,  and  systems  analysts,  we  are  often 
called  upon  to  evaluate  systems  which  service  input  traffic  and  produce  finished  products.  These 
sy  terns  are  sometimes  traditional  queues  or  networks  of  queues,  but  are  often  systems  with  queue-like 
characteristic  which  cannot  accurately  be  modeled  as  traditional  queueing  systems.  In  practice  and 
in  the  literature,  this  evaluation  is  traditionally  based  on  excercising  a  model-  of  the  service  system 
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by  subjecting  it  to  a  stream  of  input  traffic  and  estimating  or  calculating  some  expected  system 
performance  measure. 

We  feel  that  this  typical  experimental  design  is  lacking,  and  that  the  shortcomings  stem  from 
the  arbitrary  choice  of  the  distribution  of  the  input  process.  Especially  problematic  are  cases  where 
the  service  system  being  modeled  does  not  currently  exist,  where  worst-case  behavior  is  sought, 
or  where  we  wish  to  evaluate  the  system  in  situations  which  are  not  accessible  for  data  collection. 
Practical  service  system  analysis  is  interesting  only  when  the  service  system  is  in  an  environment 
where  the  workload  is  high  relative  to  the  system's  capability  to  serve.  In  all  that  follows,  we  are 
interested  in  finding  the  intensity  of  the  input  process  that  taxes  the  service  system  to  the  extremes 
of  its  capabilities,  and  using  this  intensity  as  a  measure  to  compare  systems. 

1      The  General  Service  System  Model 

The  service  systems  considered  all  have  the  following  features: 

1.  a  centralized,  controlable,  nonlattice  process  which  generates  tasks  at  a  rate  A  per  unit  time; 

2.  tasks  are  admitted  upon  generation  and  processed  by  the  system; 

3.  a  completed  task  is  ejected  from  the  system; 

4.  t ne  system  has  the  capability  to  process  as  many  as  /j.  tasks  per  unit  time. 

We  will  call  such  systems  Discrete  Event  Service  Systems  (DESSs),  see  figure  1.  We  will  allow 
the  system  to  create,  combine,  destroy,  or  absorb  tasks  -  the  DESS  need  not  be  work  conserving. 
A  DESS  must,  however,  be  mixed  or  open,  as  we  are  interested  in  how  much  work  we  can  inject 
into  the  system  without  overburdening  it.  We  may  measure  the  performance  of  the  system  using 
traditional  queueing  measures  such  as  the  number  of  tasks  resident  in  the  system  or  some  production 
cost,  or  we  may  opt  to  analyze  some  measures  which  are  particular  to  the  application  at  hand.  In 
this  work,  we  will  approach  the  work-conserving  and  nonwork-conserving  systems  seperately. 
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Figure  1:  A  simple  DESS. 


1.1      Motivating  Example 

Recently  a  mocH  was  constructed  of  all  of  the  single-channel  radio  communications  in  a  Marine 
Expeditionary  Brigade  (MEB).  A  MEB  consists  of  approximately -20,000  marines  who  perform  in 
three  subelements,  the  Ground  Combat  Element  (GCE),  the  Air  Combat  Element  (ACE),  and  the 
Combat  Servirp  S"pp°r'  Element  (CSSEV  The  single-channel  radio  network  of  the  MEB  may  consist 
of  several  thousand  radios  operating  on  several  hundred  radio  nets.  The  goal  of  the  study  was  to 
allocate  newly  purchased  antijamming  radios  to  the  different  units  in  the  MEB. 

The  radio  traffic  was  modeled  using  the  Marine  Corps'  version  of  a  structured  trafffic  model. 
There  are  sixty  classes  of  communication  task  packages,  each  representing  a  different  mission  that 
the  MEB  executes  in  battle.  Eacli  task  package,  called  a  Broad  Operational  Subtask  (BOST),  was 
composed  of  several  tasks  which  were  called  Message  Exchange  Occurrances  (MEOs).  Each  BOST 
contained  between  five  (5)  and  sixteen  (16)  MEOs,  and  each  MEO  has  a  set  of  other  MEOs  in  the 
BOST  which  must  be  completed  before  it  can  be  initiated. 

The  MEOs  are  completed  by  transmittting  a  message  from  the  specified  sender  to  each  of  the 
specified  receivers.   If  this  communication  cannot  be  accomplished  using  the  specified  doctrinal  radio 
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net,  the  sender  will  attempt  to  reach  the  receiver  using  other  nets,  relays  through  other  radios,  or 
using  messengers.  The  user  may  also  elect  to  delay  the  MEO  some  short  time  before  attempting  it 
again.  The  radios  in  the  system 

•  experience  failures; 

•  become  jammed  by  scenario-specified  jammers; 

•  hold  queues  of  MEOs  with  different  priorities; 

•  are  capable  of  moving  into  and  out  of  range  of  other  radios; 

•  can  be  used  in  voice  or  digital  mode; 

•  exit  and  join  different  radio  nets  as  they  move  around  or  attempt  to  reroute  MEOs. 

The  antijamming  radio  fails  less  frequently  and  is  harder  to  jam  than  the  old  radio.  However,  it  is 
much  more  time  consuming  to  enter  a  net  with  the  new  radio  than  with  the  old  one.  Finally,  if  an 
old  radio  tries  to  contact  a  new  radio,  the  new  radio  must  transit  into  a  less  capable  mode  in  order 
to  receive  the  MEO,  then  reenter  his  regular  net  of  new  radios. 

Our  model  features  structured  traffic  being  presented  to  a  set  of  servicing  capabilities  which  act 
semi-autonomously.  The  routing  of  the  MEOs,  sometimes  using  several  transmissions  to  accomplish  a 
single  MEO,  makes  this  model  extremely  difficult  to  analyze.  To  expedite  the  analysis  and  to  acheive 
maximum  flexibility  and  sponsor  acceptance,  an  object-oriented  simulation  model  was  constructed 
to  allow  analysts  to  build  MEB  radio  net -structures  and  test  their  capabilities  against  on  another. 

This  brings  us  to  the  search  for  the  right  rate  of  generstion  of  EOSTs  in  the  system.  This 
generation  rate  is  clearly  dependent  on  the  pace  and  nature  of  the  battle  being  experienced  by  the 
MEB.  After  some  initial  searching  through  volumes  of  data  from  the  recent  Desert  Storm  operation, 
we  found  that  the  Marine  Corps  didn't  take  time  during  their  ground  war  to  record  the  time  and  type 
of  every  BOST  they  executed!  Experience  with  the  model  showed  that  the  measured  performance 
depended  greatly  on  the  pace  of  the  presented  traffic.  When  considering  various  C3I  architectures, 
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we  found  the  ranking  of  the  radio  allocations  from  best  to  worst  changed  as  the  BOST  generation 
rate  changed.  The  sponsor  wanted  to  identify  the  best  architecture  for  the  most  intense  traffic. 

2      Work  Conserving  Systems 

We  first  study  the  behavior  of  systems  where  there  is  a  one-to-one  correspondence  between  the 
tasks  we  submit  for  processing  and  the  finished  tasks  the  service  system  produces.  Work-conserving 
queueing  models  do  not  allow 

•  tasks  to  expire  while  in  service; 

•  tasks  to  create  other  tasks  while  in  service; 

•  tasks  to  be  split  or  combined; 

•  tasks  which  never  finish  service. 

Work-conserving  queueing  system  models  are  common  in  both  the  practice  and  literature  of 
applied  probability.  In  a  typical  experiment,  we  generate  input  to  the  system  at  a  constant  rate, 
monitor  the  performance  of  the  system  either  at  fixed  intervals  or  upon  departure  from  the  system, 
and  employ  well-known  methods  of  steady-state  analysis  to  estimate  the  steady-state  average  of  the 
performance  measure. 

A  maxim  of  the  analysis  of  service  systems  is  that  the  system  will  have  stationary  long-run 
behavior  if  and  only  if  the  number  of  arriving  tasks  are,  on  average,  less  than  the  number  of  tasks 
the  system  is  capable  of  processing.  If  our  overall  system  can  work  at  a  maximum  of  /x  tasks  per 
unit  time,  we  can  input  as  many  as  //  per  unit  time  and  the  system  will  remain  stationary..  If  A  is 
our  arrival  rate  foi  the  system,  we  wish  to  manipulate  A  to  expose  /i. 

2.1      Generating  Data 

There  are  two  ways  we  can  generate  data  from  a  work-conserving  system  which  will  reveal  the 
maximum  processing  rate  in  the  system.  They  are: 
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•  input  tasks  to  the  system  at  a  rate  known  to  be  much  higher  than  the  system  can  handle; 

•  fill  the  system,  then  input  a  new  task  every  time  that  a  task  completes. 

Instead  of  choosing  a  very  high  input  rate  and  dealing  with  the  problems  of  exploding  buffer  contents 
and  a  nonrecurrent  system,  we  will  simply  close  off  the  system  and  recirculate  the  tasks  which  finish. 
This  approach  will  also  serve  as  a  good  introduction  to  the  analysis  of  systems  which  do  not  conserve 
work. 

Thus,  we  examine  a  special  kind  of  closed  queueing  network  -  one  with  a  single  loop-back  which 
all  tasks  traverse.  Let  A(r)  be  the  time-dependent  rate  of  recirculation  of  tasks  in  the  system.  So 
long  as  the  system  contains  enough  tasks  to  keep  it  working  at  capacity,  we  have  X(t)  — *  \i  as 
/  — ->  oo.  Kelly  [3]  and  Walrand  [9]  both  show  this  for  exponentially  distributed  service,  and  Disney 
and  Kiessler  [2]  make  the  extension  to  Jackson  networks.  The  result  can  be  extended  in  the  obvious 
way  by  treating  Phase-type  distributions  for  service  times,  to  produce  the  result  we  seek  (A(r)  —  // 
as  t  — ►  oo)  for  generally  distributed  service. 

Example 

We  will  demonstrate  this  method  on  the  Jackson  network  shown  in  figure  2. 

2.2      Detecting  Transition  to  Steady  State 

Let  us  simulate  the  completion  of  the  first  Ar  customers  serviced  by  the  closed  system  for  M  inde- 
pendent replications.  Let  T,  _,  be  the  j"1  time  between  recirculation  during  the  \th  replication.  Thus, 
Xj j,  i  —  1,  2, .  . . ,  M  is  a  set  of  iid  samples.  Let  Tj  =  ^t=1  TtJ/A/  be  the  average  recirculation  time 
process.  We  seek  the  index  TV*  such  that  ET,j  =  ETj  —  /j  for  all  j  >  N*.  Hence,  we  are  in  the 
setting  of  a  traditional  initial  transient  detection  problem. 
There  exist  many  ways  to  tackle  this  problem,  including 

•  cross-replication  confidence  intervals,  [10] 

•  tests  for  significant  drift; 


2     WORK  CONSERVING  SYSTEMS 


recirculation 


(5) 


(2) 


.6 


.1 


^ 


irr! 


(4) 


Figure  2:  A  Jackson  Network  designed  to  have  a  maximum  service  rate  of  0.5.  The  numbers  in 
parentheses  are  the  number  of  servers  at  eacli  station,  and  the  routing  probabilities  are  shown  on 
the  workstation  connections.  All  servers  have  unit  service  time,  and  all  buffers  are  infinite.  The 
dashed  line  shows  the  recirculation  route  added  to  force  the  system  to  serve  at  the  maximum  rate. 
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•  standardized  time  series  (STS),  [6]. 

Tti  our  pYiipripn rps    wp  Viavp  fonnH  STS  ii«pfnl    esneciallv  in  a  slichtlv  modified  version  we  have 
developed,  which  we  call  ratio  STS  (RSTS). 

2.3     Ratio  STS 

The  method  of  standardized  time  series  (STS)   [6],  produces  confidence  intervals  from  autocorrelated, 

stationary  data.     This  method  was  used  in     [7]  to  detect  the  existence  of  initialization  bias  in 

simulation  output,  and  was  sharpened  to  produce  optimal  tests  when  the  functional  form  of  the 

initialization  bias  is  known. 

Suppose  that  we  have  M  independent  samples  of  n  points  each,  with  Y,j  being  the  jth  point  in 

the  i      independent  sample.  Let 

3 

y«j  =  i"1Ey'.*  (1) 

fc=i 

for  i  =  1,2,...,  M ,  and  with  Vjo  =  0  for  each  i.  The  time  series  Si(k),  k  =  1,2, . . .,  n  is  constructed 
for  each  independent,  replication  i  as 

Y,,n  -  Y,,k     for  0  <  h  <  n 
Si(k)  =  \  (2) 

0  k  =  0,n. 

Let  a  be  the  variance  of  Y,.  If  S, (k)  is  divided  by  rr^/ri/k  and  scale  the  index  k  so  that  the  result 
resides  in  the  unit  interval  [0,  1],  the  resulting  time  series  T,(/),  0  <  /  <  1  is  known  to  approximate  a 
Browniai.  bridge  as  n  — <•  oo.  This  is  the  fundamental  result  of  [6],  and  the  theoretical  basis  of  this 
sequencial  procedure. 

Schruben  shows  that  scaling  and  summing  Ti(t), 

n 

At  =  ay/nY^Ti(kn)  (3) 

results  in  a  normal  random  variable  A,  with  variance  given  by 

im«m.)  =  ''"";.]" ]).  (4) 
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Note  that,  except  for  a  factor  of  a'1 ,  VAR(A,)  is  independent  of  the  data,  it  relies  only  on  the 
parameters  of  the  experiment.  Hence,  for  any  integer  d  <  M, 

xl  ~  ]T(      ;4;     ?  (5) 

fr[   s/VAR(Ai) 

=      ^2 -!)!>?■  (6) 

v  '  1  =  1 

The  original  STS  used  to  detect  initial  transients  used  \"d  as  a  lesl  statistic  for  stationarity  of  the 
mean  response.  If  we  form  a  ratio  of  \2d  and  \1t_d,  we  can  eliminate  the  need  to  estimate  a2, 
forming 


Fd,M-d~— „    ,  ^lf    — •  (7) 


This  test  statistic,  which  we  call  the  RSTS  test  statistic,  is  easy  to  use  in  all  of  the  applications 
where  STS  is  applied.  In  particular,  if  we  are  interested  in  determining  the  onset  of  steady  state,  we 
can  form  the  backward-moving  sequences  Aj  j,j  =  n  —  1,  n  —  2, . . . ,  1  for  each  replication  i,  where 
Ajj  is  formed  from  the  subsequence  YJ  jt,  k  =  j,j  -+■  1, . . . ,  n,  the  port  on  of  the  i  replication  between 
j  and  n.  Thus,  we  form  the  sequence  of  F-statistics 

ci~1Yd      \2 

If  we  assume  that  the  system  is  in  steady  state  when  each  of  the  A,  n  are  collected,  then  we  can  detect 
the  transition  of  the  system  into  steady  state  by  looking  at  the  first  index  Ar*  where  FdtM-d(N*) 
exceeds  the  critical  value  for  an  F  random  variable  with  identical  degrees  of  freedom.  This  method 
is  demonstrated  in  the  following  example. 

2.3.1      Example 

Continuing  with  the  work-conserving  system  example,  suppose  that  we 

•  start  the  system  with  25  tasks  enqued  at  workstation  1  at  time  0.0; 

•  simulate  N  =  500  customer  recirculations; 

•  replicate  M  —  20  times. 
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Figure  3:  Trajectory  of  the  mean  process  Tj  and  the  confidence  intervals. 


Figure  3  shows  the  trajectory  of  Tj  and  the  associated  confidence  interval  process  for  the  first 
HO  recirculations.  Clearly,  by  sample  TV*  =  80  we  have  passed  the  criteria  for  being  in  steady  state 
according  to  Welsh's  cross-replication  confidence  interval  method.  Furthermore,  we  can  see  that  any 
detectable  slope  in  the  mean  process  is  neglegible.  When  tested  for  our  20  independent  samples,  the 
drift  of  tested  to  be  insignificant  (Ho-  "o  drift  has  p-value  «  0.4). 

The  mean  time  between  recirculations  in  100,  101, .  . .,  140  was  T  —  0.56,  (confidence  interval 
(0.55330,0.56091)),  clearly  not  as  fast  as  the  \i  =  0.50  which  we  know  to  be  the  system's  capacity. 
Performing  unweighted  RSTS  in  the  first  140  samples  showed  no  transition  to  steady  state  detectable 
-  the  procedure  seemed  to  be  accurate  enough  to  discern  that  the  tiansition  had  not  yet  occured, 
see  figure  4. 
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Figure  4:  RSTS  performed  on  the  first  140  sample  recirculation  times,  using  12  numerator  degrees 
of  freedom  and  8  denominator  degrees  of  freedom.  Conclusion:  no  transition  to  steady  state  was 
evident. 
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Figure  5:  RSTS  performed  on  the  first  500  sample  recirculation  times. 


When  we  extend  the  length  of  the  runs  we  consider  to  the  full  500  samples,  we  see  that  RSTS 
was  able  to  indicate  a  strong  transition  to  steady  state  around  the  N*  =  330  sample,  see  figure 
5.  Averaging  the  samples  collected  in  350,351, . .  .,500,  we  observe  an  overall  average  of  T  =  0.501 
(confidence  interval  (0.49816,0.50389)). 

RSTS  clearly  dominated  the  other  traditional  initial  transient  methods.  In  the  case  of  the 
recirculating  jobs,  we  clearly  have  a  very  gradual  descent  to  the  steady-state  average.  The  detection 
method  is  not  important  to  our  overall  theme,  though  we  must  issue  a  general  caution:  The  choice 
of  TV*  should  be  made  very  conservatively. 
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3      Non- Work- Conserving  Systems 

In  this  section,  we  consider  the  case  where  work  is  not  conserved  in  the  system,  and  where  the 
measure  of  performance  is  very  general.  This  case,  which  is  much  more  interesting  and  applicable 
than  the  work-conserving  case,  allows  us  to  treat  cases  where  the  service  system  may  share  unusual 
characteristics  with  the  real  system,  and  where  the  measure  of  performance  of  the  system  may  be 
dictated  by  the  study  sponsor. 

3.1      Motivating  Example,  Revisited 

In  the  motivating  example,  BOSTs  were  input  to  the  system.  Depending  on  the  BOST  type,  the 
BOST  might 

•  splinter  into  several  communications  tasks,  which  may  splinter  further  at  a  later  time; 

•  require  partial  or  full  reassembly  at  different  points; 

•  expire  after  is  has  reached  a  certain  age. 

Thus,  the  system  clearly  doesn't  conserve  work. 
The  sponsor  was  interested  in  specifying 

•  a  mix  of  different  BOST  types  which  the  input  was  comprised  of; 

•  a  time  allotment  for  each  type  of  BOST,  and  a  time  when  the  BOST  expires  and  is  removed 
from  the  system; 

•  a  one-time  cost,  by  BOST  type,  assessed  when  the  time  allotment  expin  s  with  the  BOST  still 
in  process; 

All  of  these  requirements  were  imposed  because  of  the  need  to  assess  the  system's  ability  to  handle 
communications  as  diverse  as  artillery  targetting  and  mission  execution  traffic,  medical  evacuation 
requests,  situation  reports,  intelligence  traffic,  logistics  and  administrative  communications,  and 
communications  allowing  the  radio  nets  to  counter  radio  jamming. 
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For  this  application,  we  constructed  a  penalty  process  p(t),  which  superimposed  lateness  and 
expiration  penalties  from  all  of  the  traffic  in  the  system,  and  accumulated  these  penalties  in  [0,<]. 
We  compared  C3I  architectures  based  on  the  rate-of-climb  of  this  penalty  when  traffic  was  inserted 
into  the  system.  The  sponsor  really  wanted  to  know  the  answer  to  the  following  question:  "Which 
architecture  has  the  best  peak-load  performance?"  Since  the  model  was  an  idealization  of  the  real 
system,  and  since  wartime  traffic  load  data  is  not  available,  we  were  faced  with  the  problem  of 
determining  what  traffic  rate  represented  peak  loading  to  the  system. 

3.1.1      The  Workload  Ramp 

Given  that  we  do  not  know  the  service  rate  of  the  system,  we  can  try  to  produce  an  estimate 
by  modulating  the  intensity  of  the  system's  input  and  observing  the  effect  this  has  on  the  system 
performance.  One  might  attempt  this  by  stepping  through  some  reasonable  intensity  values,  or 
by  doing  some  sort  of  iterative  search.  After  considering  several  alternatives,  we  decided  that  a 
nonhomogenius  input  process  with  gradually  increasing  intensity  might  be  appropriate.  Tins  idea 
is  similar  to  testing  a  stereo  system  for  its  ability  to  play  loud  music  -  we  gradually  turn  up  the 
volume,  listening  for  the  point  where  the  music  begins  to  become  distorted. 

The  mechanics  of  generating  the  ramping  workload  process  and  calculating  the  likelihood  ratio 
for  a  generated  sample  path  are  now  presented.  The  requirement  is  to  generate  a  nonhomogeneous 
Poisson  process  with  jumppoints  x\,  £2,  ••-,%N  from  an  intensity  function  \{t)  given  by 

A(0)  +  r/     t>0 
--<  (9) 

0  t  <  0 

where  A(0)  >  0  is  the  initial  intensity  and  r  is  the  rate  of  climb  of  the  intensity  function.  In  this 
presentation  sign(i-)  is  not  specified,  and  is  a  point  of  interest  in  future  research.  Let  N(t)  be  the 
number  of  tasks  arriving  during  [0,t].  From  the  above  equation,  we  have 

a(t)     =     E[N(t)\  (10) 

=      /  [A(0)  +  rs]ds  (11) 

Jo 
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-     \(Q)t  +  ri2/2  (12) 

Thus,  if  Xi,x2,  ■  ■  ■  ,iN  is  generated  via  (9),  then  a,(x1,x2,  ■  ■  ■ ,  £jv)  is  a  Poisson  process  with  rate  1, 
[1].  The  scheme  for  generating  our  ramping  workload  process  is  given  as 

t0  =  0.0.;=  1 
while   NOT  DONE 

generate   U   ~   U[0,    l] 

Tj  <—  tj_i    -   ln(U) 


_i/       v  -A(0)+x/A(0)-'  +  2TJr 

*j  -  a      to)  =  ^7 

end   while 

Algorithm  1:  The  generation  of  the  ramping  intensity  workload  process. 

Let  G(x)  =  1  —  G(x)  be  the  complement  of  the  joint  distribution  of  i\ ,  I-j, .  .  .  ,  j\y.  Then 

Gx.lx.-^O      =      P[x,- -£,-!  >t\xt-i]  (13) 

_       e-[a(r,_1+t)-a(i«-l)]  ^4^ 

_        e-[A(0)(f1_1+()  +  r/2(f,_1+()2]-[A(0)(f,_,)  +  r/2(f,_Ir']  Qtj\ 

Thus,  the  conditional  density  function  #.?, ].?,_,  (0  is  given  by 

fc.lf.-.W      =      [A(0)  +  r(<  +  f1_1)]e-[A'0^'-/2"2+2f-'«  (16) 

=      [A(0)  +  r(;+f,_1)]f_tA(0)+r('/2+f-l)lt-  (17) 

This  last  presentation  of  <7£,|£,_,(£)  highlights  the  nature  of  the  density  function,  with  leading  con- 
stant given  as  A(0)  4-  r£,-,  the  rate  at  the  time  of  the  generation  epoch,  and  the  exponent  given  as 
-t  times  A(0)  +  r(£',_1  +  z,)/2. 

When  appropriately  calibrated,  the  ramping  workload  process  will  drive  the  DESS  into  regimes  in 
which  it  is  underutilized,  progressing  to  the  point  of  total  utilization,  and  then  becoming  saturated. 
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3.2      Detecting  the  Transition  to  Overloaded  for  a  DESS 

Let  our  system  performance  measure  for  input  intensity  A  and  time  t  have  expected  value  rp(X,t). 
Our  only  assumptions  about  V  are  that  it  is  responsive  to  changes  in  A,  it  grows  at  a  rate  similar  to 
a  degree  s  polynomial  when  A  <  ;/,  and  it  grows  faster  when  A  >  //. 

Hence,  by  estimating  the  s  +  Ist  derivative  of  the  performance  measure  (s -\-  1  because  we  would 
like  to  deal  with  a  mean-zero  sequence),  we  can  produce  a  sequence  with 

•  constant  mean  when  A  <  //, 

•  some  drift  when  A  >  //. 

Let  <i,<2i  •  ■  •  i  *n  be  tj  evenly  spaced  points  in  time,  where  A(/])  is  believed  to  be  less  than  the  service 
capacity  //  of  the  DESS,  and  A(in)  is  believed  to  be  much  greater  than  ^i.  Our  method  for  estimating 
DESS  capacity  will 

1.  replicate  the  system  performance  At  times  using  the  workload  ramp  as  the  input  process, 
collecting  data  tyl(\(tj),  tj)  for  the  ith  replication; 

2.  form  the  s  +  l"  derivative  data  vJ»|.i  +  1)(A(/j ■),  tj ■),  j  =  s  +  1,  s  +  2, . .  .,n,i  =  1,2, . .  .,n  using 
sequencial  differences; 

3.  perform  a  transition  detection  to  determine  the  point  N*  where  ty,*       (X(tj),tj)  no  longer 
have  constant  mean  0. 

We  will  propose  and  evaluate  three  methods  for  detecting  the  transition  in  the  performance 
measure: 

•  modified  A'-charts; 

•  RSTS; 

•  adaptive  regression  splines. 

The  application  of  each  of  these  three  methods  will  be  made  more  difficult  by  a  common  feature  of 
data  we  have  collected  -  although  the  mean  performance  becomes  constant  after  several  differencing 
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IT 


WORKSTATION 

P[DESTRUCTION] 

P[CREATION] 

INSERTION  STATION 

left. 

0.2 

0.0 

- 

top 

0.1 

0.4 

2 

bottom 

0.2 

0.1 

4 

right 

- 

0.3 

2 

Table  1:    Destruction  and  Creation  of  Tasks  within  the  DESS.  Tasks  that  are  destroyed  are  not. 
considered  completed. 

operations,  the  variance  still  grows.  Hence,  our  model  of  the  data  from  the  system  when  unsaturated 

will  be  Yij,i  =  1,2, ..  .,M,j  =  1,2 n,  which  are  mutually  independent,  and  where  EY  =  0  is 

constant  and  ay  is  unknown  and  assumed  to  vary. 


Example,  Continued 

We  modified  the  Jackson  network  described  in  section  2  so  that  it,  was  no  longer  work-conserving, 
and  so  that  it  exhibited  properties  similar  to  the  communications  system  described  above.  Table  1 
shows  what  can  happen  to  tasks  in  the  system  when  they  complete  service  at  each  of  the  nodes. 

For  this  system,  we  still  have  a  maximum  input  rate  of  A  =  2,  as  workstation  1  is  not  interfered 
with,  and  still  has  two  servers  serving  at  unit  speed.  The  tasks  for  the  system  are  of  three  classes. 
All  have  identical  processing  speeds  at.  the  workstations,  but  each  class  pays  a  different  price  for 
waiting  in  each  workstation  buffer.  Finally,  each  task  can  stay  in  the  DESS  cost-free  for  some  period 
of  time.  After  this  delay,  the  cost  per  unit  time  is  accumulated  based  on  the  buffer  the  task  resides 
in.  Each  task  becomes  costless  once  it  leaves  the  system.  All  of  the  data  on  task  classes  is  in  table 
2. 

The  system's  measure  of  performance  is  the  cost  accumulated  from  the  beginning  of  the  simu- 
lation. We  subjected  this  system  to  a  ramped  workload  process  which  started  with  insertion  rate 
A(0)  =    1.0  and  climbing  at  a  rate  r  =  0.00666.    Thus,  the  system  capacity  is  reached  at  time 
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CLASS 

FREE  TIME 

COST(LEFT) 

COST(TOP) 

COST(BOTTOM) 

COST(RIGHT) 

1 
2 
3 

2.0 
0.0 
1.0 

4 
1 
3 

1 

2 
2 

0 
3 

1 

1 
4 
0 

Table  2:  Free  Time  and  Holding  Costs  for  Task  Classes. 

150  =  Ar*.  We  continued  each  experiment  to  time  300. 

As  we  see  from  figure  6,  we  have  the  expected  properties  of  constant  mean  and  growing  variance 
for  the  second  derivative  of  the  accumulated  costs  when  t  <  150,  and  something  else  occurring  when 
t  >  150. 

3.2.1      Control  Charting 

The  most  straightforward  methodology  comes  from  the  field  of  quality  assurance,  [5],  and  involves 
the  use  of  a  sequencial  hypothesis  test.  We  wish  to  know  the  first  time  that  we  can  conclude  that 
it  is  no  longer  plausible  that  the  response  mean  is  EY  —  0.  Our  growing  varaince  causes  us  to  take 
one  of  two  actions: 

•  use  the  cross-replication  sample  standard  deviation  to  form  a  confidence  interval; 

•  model  the  growth  of  the  standard  deviation  and  use  the  model  when  setting  control  limits. 

Using  the  second  derivative  data  from  our  example  system,  we  can  see  that  the  former  method 
is  inclusive  because  of  constant  false  alarms  (the  lower  control  limit  makes  scleral  visits  above  the 
x-axis  during  the  trajectory),  while  using  the  modeled  standard  deviation  creates  an  envelope  which 
the  mean  stays  inside  even  when  we  know  the  system  is  exhibiting  drift,  see  figure  8. 


Let  us  return  to  the  construction  of  the  sequencial  RSTS  methods.     In  this  case,  we  have  two 
alterations  to  make: 
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Figure  6:  Trajectories  of  the  Second  Derivative  of  the  Accumulated  Cost  Samples.  The  center  line  is 
the  mean  of  M  =  20  independent  replications,  while  the  surrounding  lines  are  the  upper  and  lower 
normal  confidence  intervals. 
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Figure  7:   The  sample  standard  deviation  of  the  responses  and  a  regression  model  of  the  growing 
standard  deviation. 
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-  Figure  8:  Control  limit  detection  methods. 
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1.  the  process  must  be  reversed  to  detect  transitions  out  of  steady  state; 

2.  the  growing  variability  of  the  data  violate  the  assumptions  under  which  T,(t)  on  [0,  1]  was 
shown  to  be  a  Brownian  bridge  -  the  underpinning  of  the  method. 

Let  the  sequences  Aij,j  =  1,2,  .  .  .,j  for  each  replication  i,  where  Al<3  is  formed  from  the  sub- 
sequence Y,  k  k  =  1,2,..., j.  Thus,  we  are  moving  through  the  output  sequence  in  the  forward 
direction  (opposite  the  usual  STS  for  initialization  bias). 

The  growing  variability  of  the  output  must  be  attacked  directly.  The  model  which  fits  the  output 
with  A  <  //  is  a  mean-zero  model  with  linearly  increasing  standard  deviation.  The  expected  number 
of  input  tasks,  a(t),  is  also  quadratically  increasing  and  is  intimately  linked  to  the  growth  of  the 
preformance  measure  and  its  variability.  If  we  collect  samples  of  the  cost  function  0(A,<)  such 
that  there  are  a  constant  expected  number  of  input  tasks  per  sample,  we  can  avoid  the  problem  of 
growing  variability.  Our  data  will  still  be  formed  by  taking  sequencial  differences,  but  the  intervals 
of  sampling  will  contract.  Figure  9  shows  empirically  that  sampling  this  way  produces  data  with 
constant  standard  deviation  in  our  example.  Proving  that  this  technique  works  in  all  cases  is 
impossible  because  of  the  breadth  of  our  generality  here. 

Let  us  establish  the  time  interval  sequence  ti^2,  ■  ■  ■  ,tn  such  that  we  expect  to  inject  exactly 
C  tasks  into  the  system  between  /,  and  £,-+i,»  =  1,2,...,?!  —  1.  Hence,  given  tt  we  can  compute 
£j+i  =  t,  -f  At  using 

a(ti  +  At)  -  afc)  =  C.  (18) 

Using  12,  we  produce  the  equation 

X(0)At +  ^-(2ttAt  +  A/2)-C  =  0,  (19) 

leading  to 


±t  _  -(rt,  +  A(0))  +  jrHj  +  2r/,A(0)  +  A(0)?  +  2KJ 

r 

If  we  sample  the  DESS  cost  function,  starting  with  some  initial  point  t\  and  with  some  arbitrary 
value  C,  we  produce  a  sample  with  constant,  standard  deviation.  Taking  the  sequencial  differences 
as  above,  we  produce  a  sequence  which  has  the  properties  required  to  perform  RSTS. 
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Figure   9:    Sample  standard   deviations  from   data   collected   using   the   modified   sampling  points 
method. 
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Figure  10:    Trajectory  for  RSTS  for  the  modified  sampling  points  data  collection.    Result:    clear 
transition  at  time  123.0 

In  our  example  system,  we  used  the  modified  sampling  points  to  produce  the  F-statistic  trajectory 
shown  in  figure  10.  The  RSTS  method  gave  a  clear  signal  that  the  system  lost  steady-state  at  time 
123.0,  where  \(t)  =  1.82  and  p  —  0.910.  We  also  tried  RSTS  with  evenly  spaced  points,  and  acheived 
approximately  the  same  result.  From  this  experience,  we  conclude  that  straightforward  RSTS  is 
fairly  robust  with  respect  to  fluctuations  cr  growth  in  variability  when  the  changes  are  coordinated 
as  they  are  in  this  experiment.  Analytical  investigation  has  shown  that  these  fluctuations  do  not 
cancel  in  the  RSTS  test  statistic. 
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3.2.3      Adaptive  Regression  Splines 

Adaptive  regression  splines,  especially  multivariate  regression  splines,  are  the  focus  of  intense  basic 
research,  see  [8].  In  our  application,  the  regression  spline  required  is  especially  simple,  a*  the  model 
is  of  a  single  independent  variable,  and  we  seek  a  single  knot  in  the  regression  spline  at  the  point 
where  the  zero-mean  model  departs  from  the  data.  Using  the  methods  in  Larson  [4],  we  can  derive 
the  location  of  this  single  knot  analytically. 
The  model  used  is  stated  as 

Y(t)=\  (21) 

(    0o  +  Mt-  N*)  +  U     t  >  N', 

where  N"  is  still  our  point  where  the  regression  model  changes.  If  we  further  prescribe  that 

0o  =  0,  (22) 

ft  =  0,  (23) 

we  will  be  fitting  the  model  with  zero  mean  to  the  left  of  the  single  knot.  Let 

ssel  =  J2  y2(0;  (24) 

t<N' 

SSER=   J2(Y(t)-p2(t-N*))2\  (25) 

t>N' 

SSE  =  SSEL  +  SSER .  (26) 

The  method  involves  the  optimal  location  of  A^*  to  minimize  SSE.  The  procedure  provided  in  [4] 
can  be  simplified  for  our  application  into  a  single-pass  examination  of  the  means  of  the  data,  but  is 
valid  only  when  a  linear  model  is  appropriate  for  the  data  to  the  right  of  N* .  This  method  is  valid 
when  there  is  growing  variability,  as  seen  in  our  application. 

In  our  example,  we  fit  the  model  in  (21)  with  three  parameters,  then  restricted  it  to  the  mean- 
zero  model  prescribed  in  (22-23).  Figure  11  shows  the  two  regression  models.  The  three  parameter 
model  located  a  knot  at  time  Ar*  =  136.25,  when  the  input  rate  is  1.91  and  the  traffic  intensity  is 
0.95.  With" the  one  parameter  restriction,  the  adaptive  spline  located  the  optimal  knot  at  128.0.  We 
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Figure  1 1 :  Adaptive  regression  spline  with  a  single  knot,  three  parameter  and  one  parameter  models. 

should  note  that  the  linear  model  does  fit  fairl>  well  to  the  right  of  A7*,  as  one  might  expect  of  any 
queuing  application  with  less  capacity  than  required. 

4      Conclusion 


In  this  work,  we  have  described  the  various  problems  which  arise  when  we  are  attempting  to  de- 
termine the  service  capacity  of  a  black  box  service  system.  We  divided  this  investigation  into  two 
distinct  parts,  one  dealing  with  simple  queueing  systems  which  conserve  work.    In  this  case,  we 
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showed  the  advantages  of  closing  the  system  so  that  output  from  the  system  was  recirculated,  as  the 
time  between  recirculations  converges  to  the  system's  service  rate.   We  investigated  ways  to  detect 
this  convergence,  and  showed  how  difficult  this  is  in  a  simple  Jackson  network  example. 
The  second  part  of  the  exploration  dealt  with  systems  which 

•  do  not  conserve  work; 

•  have  their  performance  measured  using  a  general  holding  cost  mechanism  or  some  other  per- 
formance measure. 

In  this  case,  we  showed  that  if  we  modulated  the  input  process  using  a  ramped-intensity  workload 
process,  we  could  drive  the  system  from  underutilized  to  saturated.  We  explored  three  methods 
which  unveil  the  point  where  this  transition  takes  place,  and  demonstrated  each  on  an  example. 

The  wider  significance  of  this  work  is  the  beginning  of  an  exploration  for  empirical  methods 
for  determining  the  capacity  of  a  service  system.  This  exploration  is  done  not  by  representing  the 
system  using  a  queuing  model  which  we  know  how  to  analyze  a  priori,  but  by  using  a  realistic  model 
of  the  system  and  measuring  its  performance  in  terms  the  user  has" in  mind. 
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